High-Resolution Simulations of Convection Preceding Ignition in 
Type la Supernovae Using Adaptive Mesh Refinement 

A. Nonaka^ A. J. Aspden^'^ M. Zingale^ A. S. Almgren^ J. B. Bell^ S. E. Woosley^ 

ABSTRACT 

We extend our previous three-dimensional, full-star simulations of the final 
hours of convection preceding ignition in Type la supernovae to higher resolu- 
tion using the adaptive mesh refinement capability of our low Mach number code, 
MAESTRO. We report the statistics of the ignition of the first fiame at an effec- 
tive 4.34 km resolution, and general fiow field properties at an effective 2.17 km 
resolution. We find that off-center ignition is likely, with radius of 50 km most 
favored and a likely range of 40 to 75 km. This is consistent with our previous 
coarser (8.68 km resolution) simulations, implying that we have achieved suffi- 
cient resolution in our determination of likely ignition radii. The dynamics of the 
last few hot spots preceding ignition suggest that a multiple ignition scenario is 
not likely. With improved resolution, we can more clearly see the general fiow 
pattern in the convective region, characterized by a strong outward plume with 
a lower speed recirculation. We show that the convective core is turbulent with 
a Kolmogorov spectrum and has a lower turbulent intensity and larger integral 
length scale than previously thought (on the order of 16 km s~^ and 200 km, 
respectively), and we discuss the potential consequences for the first fiames. 
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1. Introduction 

For the Chandrasekhar mass white dwarf (single-degenerate) progenitor model of Type 
la supernovae (SNe la), the location of the first flames greatly affects the outcome of the 
explosion (see for example Niemeyer et al. 1996; Plewa et al. 2004; Livne et al. 2005; Garcia- 
Senz & Bravo 2005). The convective state leading up to ignition is highly nonlinear and the 
ignition results from a hot temperature perturbation near the center of the white dwarf. Once 
the temperature exceeds ~8 x 10^ K, a hot spot burns faster than it can cool via expansion 
(Nomoto et al. 1984), igniting a flame. In earlier studies on white dwarf convection in SNe 
la (Zingale et al. 2009, henceforth Z09; Zingale et al. 2011, henceforth Zll), we performed 
three-dimensional, full-star simulations of the flnal ^3 hours of convection in a white dwarf 
leading up to the ignition of the flrst flames. We followed the nonlinear rise in the temperature 
approaching ignition and showed that the ignition is likely to take place off-center (50 km 
most favored, with a likely range of 40 to 75 km, and an outer limit of 100 km) in an outward 
flowing parcel of fluid. Our results differed from the two-dimensional wedge simulation of 
Hoflich & Stein (2002) which argued that the ignition is closer to the center (^ 30 km), and 
is driven by inflow compression. 

It is important to understand how robust our results for the likely ignition radius are 
to resolution. With the recently implemented adaptive mesh reflnement (AMR) capability 
in our low Mach number code, MAESTRO (Nonaka et al. 2010, henceforth NIO), we are 
now able to study the flnal minutes of convection up to ignition at unprecedented resolution. 
We are also interested in the likelihood of multiple ignition points. Detailed visualizations 
of the evolution of the last few hot spots preceding ignition will be used to examine this 
scenario. Previous studies with an anelastic code showed that a dipole flow dominates the 
flow (Kuhlen et al. 2006, also seen in the follow-up studies shown in Woosley et al. 2007; 
Ma 2011). Our results for non-rotating white dwarfs also saw this feature. Here we examine 
this structure at higher resolution. 

Higher resolution is also important for resolving the turbulence and capturing the tur- 
bulent cascade. Simulations have shown that the flame needs to accelerate considerably 
beyond its laminar value for the resulting energetics to match observations. The primary 
mechanism for this acceleration is thought to be instabilities and the interaction with tur- 
bulence (Mueller & Arnett 1986; Livne 1993; Khokhlov 1995; Niemeyer & Hillebrandt 1995; 
Niemeyer & Woosley 1997). A popular view is that the flame interacts with turbulence 
generated by the flame itself via instabilities. The vast majority of simulations to date have 
only considered this flame-generated turbulence during the explosion phase. Aspden et al. 
(2011) suggested that turbulent entrainment was the dominant mechanism for enhancing 
the burning rate, and that the local flame speed, whether laminar or turbulent, was largely 
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unimportant. As the flame encounters lower densities and broadens, the turbulence may 
be able to directly affect the flame structure (at this point, the flame is said to be in the 
"distributed burning regime"). The altering of the flame by turbulence has been the subject 
of many studies, both semi-analytic and one-dimensional calculations with a model for tur- 
bulence (Lisewski et al. 2000; Pan et al. 2008; Woosley et al. 2009) and multi-dimensional 
numerical simulations (Ropke et al. 2004; Aspden et al. 2008a, 2010, 2011). If the conditions 
are right, the flame may transition to a detonation in this regime (Khokhlov et al. 1997; 
Niemeyer & Woosley 1997; Woosley et al. 2009, 2011). 

What are not well known are the characteristics of the turbulence that already exists 
at ignition from the centuries-long convective period. The very flrst flame (s) that ignite will 
form flame "bubbles" that buoyantly rise away from the center as they burn outward. These 
bubbles will deform due to shear instabilities and interact with the pre-existing turbulence 
and wrinkle (Garcia-Senz & Woosley 1995; Bychkov & Liberman 1995; lapichino et al. 2006; 
Zingale & Dursi 2007; lapichino & Lesaffre 2010; Aspden et al. 2011). If the turbulence is 
strong enough, it could potentially disrupt the flames or even quench them. Additionally, 
the initial convective velocity fleld has been shown to introduce large asymmetries in the 
burning (Livne et al. 2005). 

In this paper, we expand upon our previous studies of the final hours of convection 
leading up to the ignition of the first fiames in Type la supernovae. In Z09, we used 13.2 km 
resolution; in Zll, we used 8.68 km resolution. Here, we use the AMR capability of MAE- 
STRO to compute ignition statistics at 4.34 km resolution, and general fiow field properties 
at 2.17 km resolution. 

This paper is organized as follows. In Section 2 we give an overview of the MAESTRO 
algorithm, including our latest improvements for both regridding and adding an additional 
level of refinement to a simulation in progress. In Section 3, we describe our new high- 
resolution simulations. We examine the ignition statistics and compare them to our previous 
results in order to show that we have achieved suflicient resolution in our determination of 
likely ignition radii. We determine the likelihood of multiple ignition points by examining 
the dynamics of the last few hot spots leading up to ignition. We provide visualizations of 
the convective flow fleld to gain a better understanding of the flow structure. We include a 
detailed analysis of the turbulent nature of the flow fleld, and discuss the implications for 
the flrst flames. Finally, in Section 4, we summarize and conclude. 
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2. Numerical Methodology 

MAESTRO is a finite- volume, AMR hydrodynamics code for low Mach number astro- 
physical flows. In our low Mach number formulation, sound waves have been analytically 
removed, allowing for a time step based on the fluid velocity CFL constraint rather than 
the sound speed CFL constraint, while retaining compressibility effects due to background 
stratiflcation, reaction heating, and compositional changes. The algorithm is described in 
full detail in NIO. We note that the low Mach number equations do not enforce that the 
Mach number remain small; rather, if the dynamics of the flow are such the Mach number 
does remain small, then these equations are valid approximations for the evolution of the 
flow. Thus, MAESTRO is not suitable for post-ignition calculations, where we expect the 
Mach number to quickly approach or exceed unity. Also, our low Mach number approach 
assumes that the background state is spherical; thus, any deformation due to rotation is not 
accounted for in the background state. 

We now summarize the algorithm, and then describe the new procedures for dynamically 
changing the grids as well as adding an additional level of reflnement to a simulation in 
progress. For the simulations in this paper, we begin with data from Zll, in which we 
computed the last ^3 hours of convection in a non-rotating white dwarf up to the point of 
ignition using 8.68 km resolution (576^ computational cells; the problem domain is 5000 km 
on a side) and no AMR. We expand upon this simulation by adding a level of reflnement a few 
minutes before ignition and examining the ignition statistics for a 4.34 km (1152^ effective 
grid cells) resolution simulation. Next, we will add an additional level of AMR to examine the 
turbulent flow fleld in a 2.17 km (2304^ effective grid cells) resolution simulation. Computer 
allocation limits prevent us from running 2.17 km resolution simulations to ignition, even 
with the efficiency gains provided by AMR. 

2.1. MAESTRO Overview 

MAESTRO is based on the BoxLib software framework (Rendleman et al. 2000), which 
provides infrastructure for block-structured AMR applications, and includes linear solvers 
that scale to 100,000 cores on the current generation of supercomputers (see Almgren et al. 
2010 for details). We use a finite- volume approach, where each computational cell stores 
the average value of a state variable in that cell. The domain is decomposed into a nested 
hierarchy of logically rectangular Cartesian grids with computational cell width Ax^ in each 
direction (the grids at the coarsest level are associated with level index ^ = 1, the first level 
of refinement with ^ = 2, etc.), and a refinement ratio of two in each spatial direction. We 
solve a system of coupled PDEs containing advection and reaction terms constrained by an 
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equation of state that takes the form of a divergence constraint on the velocity field. This 
constraint is enforced using a projection method, which requires linear solvers to solve a 
variable-coefficient Poisson equation. 

One feature that makes MAESTRO different from standard AMR hydrodynamics codes 
is the presence of base state variables, which are functions of radius and time, (r, t), as 
opposed to Cartesian grid quantities which are functions of all spatial dimensions and time, 
(cc, t). We represent base state variables using a one-dimensional, time-dependent array. The 
base state array has a constant grid spacing, Ar = Ax^™^^/5, where ^^ax is the ffnest level 
in the simulation, and due to the spherical nature of our problem, does not directly align 
with the Cartesian grid. Figure 1 shows a depiction of the Cartesian grid, one-dimensional 
radial array, and a graphical representation of how they relate to each other. Some base state 
variables are cell-centered, and others are defined on edges. Each of the base state variables 
is computed directly from other base state variables and/or Cartesian grid variables. The 
base state density obeys an evolution equation within each time step (described below). We 
require frequent mapping from the base state to the Cartesian grid, and vice versa. In NIO, 
we describe how we interpolate a base state quantity onto the Cartesian grid, as well as a 
"lateral average" procedure that determines the average value of a Cartesian grid quantity 
at a particular radius and maps that value onto the radial array. 

In the following overview, all variables are assumed to live on the Cartesian grid, unless 
noted otherwise. MAESTRO solves the equations of reacting fiow constrained by an equation 
of state in the form of a divergence constraint. The species are evolved according to 

= -V • (pX.XJ) + poj,, (1) 

where p is the density, is the mass fraction of species fc, U is the velocity field, and ujk is 
the creation rate of species k due to reactions. We note that the density can be determined 
at any time using 

p = J2ipXk), (2) 

k 

and thus density does not have to be explicitly evolved in time. 

We define a base state density, po{r^t)^ that represents the average value of density at 
a particular radius. The base state density has its own evolution equation, as described 
below. The base state (thermodynamic) pressure, ]9o(r, t), is computed using the condition 
of hydrostatic equilibrium, 

Vi?o = ^ = -Pog^ (3) 

where the gravity, g{r^t) is computed by integrating po assuming piecewise-constant profiles 
of po within each radial cell. 
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In general, given p and X^^ we could derive the temperature from the specific enthalpy, 
/i, evolved as 

^ = -v. (pro) + ^ + (4) 

where i/nuc is the energy generation rate from reactions. In practice, we adopt the prescrip- 
tion used in Z09, Zll and derive the temperature from p, X/^, and po, effectively decoupling 
the enthalpy from the problem. In the future, we will seek ways to evolve the enthalpy in a 
manner that minimizes the drift from the equation of state. 

The velocity field is evolved according to 

|^ = -U.VU-iv.-^^je,., (5) 

ot p p 

where tt is the perturbational pressure, i.e., the local deviation of the total pressure from po, 
and Cr is the unit vector in the outward radial direction. The evolution of the thermodynamic 
variables (p, X^, and po) are constrained by the equation of state, which we represent as a 
divergence constraint on the velocity field, 

V.(AU)^A(s-=i-f). (6) 

Here, /3o{r^ t) is a base state quantity that captures the expansion/contraction of a fluid parcel 
as it changes altitude, and is a local source term that captures the compressibility effects 
due to reactions and compositional changes. The quantity Ti{r^t) is a base state variable 
representing the average at constant radius of Fi = d\ogp/d\ogp\s^ where s is the entropy. 
A full derivation of this constraint, the form of /3o and and the numerical projection can 
be found in Almgren et al. (2006a, b); Almgren et al. (2008). 



The evolution equation for po is 

5po _ d (poWq) dr]p 



dt dr dr 

where w^{r^t) is the base state expansion velocity, which accounts for the expansion of 
the atmosphere due to large-scale heating. We compute this term by integrating a one- 
dimensional version of the divergence constraint (Eq. [6]). The quantity /^^(r, t) is a base state 
quantity that accounts for changes in background stratification due to large-scale convection 
(see Almgren et al. 2008 and NIO). 

The velocity field can be decomposed into the base state velocity and a local velocity, 
U(a3,t), that governs the local dynamics, 

\]{x, t) = wo{r, t)er + U(a^, t). (8) 
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We follow the approach in NIO, where we compute the evolution of these terms separately, 
and thus evolve U subject to a perturbational form of Equations (5) and (6). 

We note that the base state quantities po^Po^ /^o^ Ti, and r]p are stored on radial cell- 
centers, whereas Wq is stored at radial edges. 

To summarize, we advance Equations (1), (5), and (7), subject to Equations (2), (3) and 
(6). We use a second-order predictor-corrector approach in which we discretize the advection 
terms using an unsplit Godunov method, compute the effect of reactions on a cell-by-cell basis 
using the VODE stiff ODE package, and couple these processes using Strang splitting. We 
enforce the divergence constraint on velocity using a projection method, which uses multigrid 
to solve a variable-coefficient Poisson equation for the update for the perturbational pressure, 

TT. 



2.2. Regridding and Adding a Level of Refinement 

Regridding is the process of redefining the AMR grid structure based on user-specified 
refinement criteria. The regridding algorithm also uses interpolation stencils to initialize 
data on newly created refined grids from underlying coarse data. Here we have improved 
the regridding algorithm described in NIO and have also implemented a new algorithm 
for adding an additional level of refinement to a simulation in progress. Our approach 
to AMR uses a nested hierarchy of logically rectangular grids with successively finer grids 
at higher levels. This is based on the strategy introduced for gas dynamics by Berger & 
Colella (1989), extended to the incompressible Navier-Stokes equations by Almgren et al. 
(1998), and extended to low Mach number reacting flows by Pember et al. (1998) and Day 
& Bell (2000). We refer the reader to these works for more details. The complication in 
applying these methods to MAESTRO is the presence of the time-dependent base state 
variables. We refer the reader to NIO for the MAESTRO-specific implementation including 
details on creating and managing the grid hierarchy, communication between levels, and the 
implementation of AMR with time-dependent base state variables. 

We note an error in the Cartesian grid regridding procedure as described in Section 5.1 
of NIO. For problems with three or more total AMR levels, we require that each grid at 
level ^ + 1 be a distance of at least four (not two as previously reported) level i cells from 
the boundary between level £ and level £ — 1 grids; this allows us to always fill "ghost cells" 
at level £ + 1 from the level £ data (or the physical boundary conditions, if appropriate). 

The major change regarding the regridding of the Cartesian grid data is in the way we 
interpolate coarse data to fill newly created fine grids. Our piecewise-linear interpolation 
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algorithm applied to pX^ causes an artificial buoyancy term to appear in the momentum 
equation, leading to spurious velocities emanating from the coarse-fine interface. The basic 
idea of the improved algorithm is to interpolate p' and Xj. separately, rather than pXj.^ to 
initialize data on the newly created refined grids. 

The variables on the Cartesian grid that we need to interpolate are U, p, pX^, Vtt, and 
S. The base state does not change structure, but we still need to recompute po, J^o, /^o, and 
Fi to be consistent with the data on the Cartesian grid. We keep the original values of w^) 
and Tjp. Here are the steps for regridding: 

1. Starting with level 1 and user-defined refinement criteria, tag all level 1 Cartesian cells 
that satisfy the criteria for refinement. Create the level 2 grids, and initialize the 
level 2 data by copying from the previous grid structure where possible. Otherwise, 
use piecewise linear interpolation from underlying coarse cells to initialize any newly 
created refined regions, including ghost cells. Continue to add additional levels of 
refinement in this way until all data on the grids at level ^^ax are filled in. There is a 
slight modification to the interpolation procedure for pXj.^ where we first interpolate 
p' = p — Pq and Xk = {pX^)/ p to newly refined regions and then construct pX^ = 
{p' + PQ)Xk. 

2. Recompute po by calling the lateral average routine, then use Equation (3) to compute 
Po- 

3. Recompute T and Fi on each Cartesian cell using the equation of state. Recompute 
Fi by calling the lateral average routine. Then, recompute /3o as described in NIO. 

4. Compute a new appropriate At. 

The procedure for adding an additional level of refinement to a simulation in progress 
is largely based on the standard regridding procedure, except that now the base state array 
will have twice as many cells since Ar is based on the resolution of the finest Cartesian grid, 
i.e., Ar = Ax^™/5. To add one additional level of refinement to a simulation in progress: 

1. Perform Step 1 in the regridding procedure defined above, except allow for an additional 
level of refinement. 

2. Define a new base state array with twice the resolution, i.e., Ar = Ax^™/5. Call the 
lateral average routine to compute po and use Equation (3) to recompute po. 



3. Recompute T and Fi on each Cartesian cell using the equation of state. Call the lateral 
average routine to compute Fi. Then, recompute /3o as described in NIO. 
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4. The base state variable Wq is edge-centered. We compute Wq on the finer base state 
array using direct injection from the previous coarser base state array on ahgning edges, 
and piecewise-hnear interpolation on non-aligning edges. 

5. The base state variable ijp is cell-centered. We interpolate rjp onto the finer base state 
array using piecewise-linear interpolation from the previous coarser base state array. ^ 

6. Compute a new appropriate At. 



3. Results 



We now focus on one particular simulation performed in Z 11, in which we computed the 
last ^3 hours of convection up to the point of ignition for a non-rotating white dwarf using 
8.68 km resolution (576^ computational cells) and no AMR. As before, we define ignition 
as the time when the maximum temperature exceeds 8 x 10^ K. Here is a summary of our 
results from that simulation: 



• The plots of peak temperature, peak radial velocity, and peak Mach number as a 
function of time each exhibited a gradual, non-linear rise until the peak temperature 
exceeded ^7 x 10^ K. Then, the rise in each field became much steeper, with ignition 
following shortly afterwards. 

• The first cell to ignite was 25.7 km off-center, and had an outward radial velocity of 
5.1 km s~^. 

• For the last ^3 minutes preceding ignition, the average radius of the hottest cell was 
52.3 km with a standard deviation of 25.5 km. 

• Histograms of the radius of the hottest cell during the final ^3 minute preceding 
ignition averaged over small time intervals indicated that 

— The favored ignition radius was 50 km, with a likely range of 40 km to 75 km, 
and an outer limit of 100 km. 



^In practice, we store r]p on radial cell centers and edges as separate arrays. We interpolate the radial 
cell-centered and edge-based arrays onto the finer base state arrays separately, rather than simply interpolate 
the radial cell-centered array onto the finer base state array, and then arithmetically average to get the radial 
edge-centered array. In the future, we will run in the latter mode for simplicity, noting that the effects of 
this change are very minor, and that both methods are second-order. 
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— Nearly all of the hot spots had an outward radial velocity. 

— These two results were consistent within any smaller time window within the final 
^3 minutes. 

• Visualizations of the convective flow fleld showed a dipole structure in the interior 
convectively unstable core, and a sharp boundary between the interior and the stably 
stratifled outer portion of the star. 

In this section, we examine the robustness of the ignition results at higher resolution. 
Then, we use new visualization techniques to follow the last few hot spots preceding ignition 
to determine the likelihood of multiple ignition points. We also visualize the overall convec- 
tive flow fleld to show the detailed flne-scale structure, as well as a more coherent picture of 
the large-scale features. Finally, we analyze the turbulent spectrum to quantify the extent to 
which we have resolved the turbulent cascade, and discuss the eflFect that turbulence could 
have on the flrst few flames. 

We note that we do not consider a high-resolution rotating white dwarf at this time. 
A 576^ rotating simulation developed high velocities on the surface of the star at the poles, 
likely due to the deformation of the star. In our lower resolution rotating runs in Zll, we 
saw a similar feature, but the velocities did not become large enough to restrict our time 
step as they do for the higher-resolution runs. A potential future solution to this would be to 
reformulate the base state in MAESTRO to deal with equipotentials instead of a spherical 
radius. 



3.1. Problem Setup 

The 8.68 km resolution simulation in Zll followed the last ^10500 s preceding igni- 
tion. The simulation required ^6 million CPU hours on the Jaguarpf Cray XT5 at OLCF. 
Assuming perfect scaling and no AMR overhead, a 4.34 km simulation from t = would 
require a factor of ^4 more computational resources (since the time step is a factor of two 
smaller, and with our tagging criteria we have approximately the same number of cells at 
levels 1 and 2, so there are twice as many total grid cells). Due to computer allocation 
limits, running 4.34 km resolution from t = is infeasible, so instead we add an additional 
level of reflnement to an 8.68 km simulation at a time corresponding to ^250 s preceding 
ignition. The edge of the star lies where po ^ 1 x 10^ g cm~^, corresponding to a radius of 
r ^ 1890 km. We reflne all level 1 cells where p > 5 x 10'' g cm~^ (r ^ 1225 km), which more 
than encompasses the convective region (the convective region boundary lies approximately 
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where po ^ 1-26 x 10^ g cm~^, with r ^ 1030 km). This new simulation has 4.34 km resolu- 
tion (1152^ effective grid cells). We note that since this problem is highly non-linear, we do 
not expect the ignition to occur at exactly the same time. In fact, the 4.34 km simulation 
takes ^350 s to ignite. Approximately 100 s into the 1152^ simulation, we add another level 
of refinement, tagging all level 2 cells where p > 1 x 10^ g cm~^ (r ^ 1080 km). This sec- 
ond new simulation has 2.17 km resolution (2304^ effective grid cells). We run the 2.17 km 
simulation for ^80 s, and not to ignition (again due to computer allocation limits). 

The resulting three-level grid structure is shown in Figure 2. The grids adaptively 
change as the simulation progresses, but since the overall base state density profile of the 
star is relatively constant (as shown in Zll), the grids do not change significantly over time. 
Some specific details concerning this grid structure are as follows. 

• The red grids are at 8.68 km resolution. There are 1728 red grids, each of which has 
48^ grid cells (^191 million grid cells total). 

• The green grids are at 4.34 km resolution. There are 1736 green grids of varying size, 
with a maximum of 48^ cells per grid. (^141 million grid cells total). 

• The blue grids are at 2.17 km resolution. There are 3646 blue grids of varying size, 
with a maximum of 64^ cells per grid. (^654 million grid cells total). 

By contrast, a simulation without AMR at 2.17 km resolution would contain 2304^ = 12.2 
billion grid cells, or a factor of ^12 more grid cells than the AMR simulation. 

For the simulations in this paper, we use the recently implemented hierarchical approach 
to parallelism described in Almgren et al. (2010). We use a coarse-grain parallelization 
strategy to distribute grids to nodes, where the nodes communicate with each other using 
MPI. We also use a fine-grain parallelization strategy in the physics-based modules and the 
linear solvers, in which we use OpenMP to spawn a thread on each core on a node. Each 
thread operates on a portion of the associated grid. Grids at each level of refinement are 
distributed independently. This approach allows for MAESTRO (in particular the linear 
solvers) to scale to ^100, 000 cores. All simulations were performed on the Jaguarpf Cray 
XT5 at OLCF using 1,728 MPI processes with 6 threads per MPI process (10,368 total 
cores) . 
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3.2. General Behavior 

We begin by reproducing some of the diagnostics used in Zll using data from the 
4.34 km and 2.17 km simulations. In Figure 3, we plot the peak temperature as a function of 
time for the 4.34 km and 2.17 km resolution simulations, and also include original 8.68 km 
resolution data for comparison. First, we see that over the last few minutes, the temperature 
profiles have the same general trend. The peak temperature value steadily grows with time, 
including fiuctuations of several percent. Once the star ignites, the peak temperature rapidly 
increases beyond 8 x 10^ K. We consider the local temperature peaks preceding ignition to 
be "failed" ignition points, i.e., hot bubbles that are not quite hot enough to cause the 
temperature to run away. The ignition time for the 8.68 km and 4.34 km simulations differ 
by ^100 s. Due to the highly non-linear nature of this problem, this result is not particularly 
surprising. At the beginning of the 4.34 km simulation, we notice that the peak temperature 
curves track each other very well for the first ^80 s (from time range 10200-10280 s) before 
the curves begin to show different behavior. This is not particularly surprising either since 
the 4.34 km solution begins as an interpolated imprint of the 8.68 km simulation. After 
^80 s, we say that the peak temperature in the 4.34 km simulation has decorrelated from 
the 8.68 km simulation, and we expect the statistical properties of the hot spots near the 
center of the star to be consistent with an independent 4.34 km simulation initialized at time 
t = 0. We still expect the general convective fiow field to look qualitatively similar for a 
longer period of time. The 2.17 km simulation shows similar behavior; after initializing the 
simulation from the 4.34 km data, it takes ^40 seconds for the peak temperature curves to 
decorrelate. 

We would like to comment on the time step used in these simulations. Using an ad- 
vective CFL number of 0.7, the average time steps over the time range 10300-10380 s are 
approximately 0.027 s (for the 8.68 km simulation), 0.016 s (4.34 km), and 0.010 s (2.17 km). 
The time steps do not quite change by a factor of two with refinement since the peak veloci- 
ties do not necessarily lie in the refined convective region. We also want to comment on the 
efficiency of the low Mach number formulation as compared to an explicit, fully compressible 
approach. In Nonaka et al. (2011), we showed that the 8.68 km simulation took a time step 
of a factor of ^70 larger than a compressible code, yet a low Mach number time step takes 
approximately 2.5 times as long given the same computational resources, yielding an overall 
efficiency increase of ^28 over a compressible code. This comparison is especially meaningful 
because we compared to the CASTRO (Almgren et al. 2010) compressible code, which is 
based on the same BoxLib framework as MAESTRO, and uses the same unsplit Godunov 
advection formulism, same equation of state, and same reaction network ODE solver. 

Next, in Figures 4 and 5 we plot the peak Mach number and peak radial velocity as 
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a function of time. We see the same general behavior as in the 8.68 km simulation, where 
the peak Mach number and radial velocity remain relatively constant until the final seconds 
preceding ignition, where the values rapidly increase. 

The next quantities of interest are the radius of the first ignited cell and its corresponding 
outward radial velocity. In the 4.34 km simulation, the radius of the ignited cell is 41.3 km, 
with outward radial velocity of 9.5 km s~^. We compare these values to those reported in 
Table 1 of Zll; the 8.68 km simulation had an ignition cell radius of 25.7 km with outward 
radial velocity of 5.1 km s~^. 

We would also like to examine the ignition radius and radial velocity if we were to define 
ignition as 1.1 x 10^ K rather than 8 x 10^ K. However, by advancing the solution using 
our advective CFL condition, the simulation quickly becomes unphysical. Specifically, if we 
continue to let the white dwarf evolve past 8 x 10^ K, over the next ^0.5 s (^60 time steps), 
the peak temperature steadily climbs to ^9 x 10^ K while the peak Mach number holds 
steady at ^0.1. Then, over the next few time steps, the temperature unphysically spikes to 
^8 X 10^ K, with the peak Mach number quickly climbing to well over 1000. Our low Mach 
number model has obviously broken down, so these results are not physical. To remedy 
this situation, and to advance our simulation to 1.1 x 10^ K, we apply a heuristic time step 
limiter, which attempts to reduce the time step so the peak temperature does not increase 
by more than ~1% each step. We limit the time step using 



At = min 



100 

J- jy^ax max 



(9) 



where A^cfl is the time step computed using our standard advective CFL condition, T^^^x 
the maximum temperature in the domain at the current time step, and T^~^ is the maximum 
temperature in the domain from the previous time step. In doing this, we find that that we 
reach 1.1 x 10^ K at 0.57 s after 8 x 10^ K, the ignition point has advected to a larger radius 
(48.9 km), and the ignition point radial velocity has increased to = 14.0 km s~^. These 
results are not surprising, given the ignition conditions at 8 x 10^ K reported above. 

In Zll we studied the time history of the hottest cell over the last few minutes. We 
gathered statistics to help us in our determination of likely ignition radii, and repeat the 
same diagnostics here. In Figure 6, we show the radius of the hottest cell as a function of 
time for the final seconds preceding ignition for the 4.34 km simulation. In Table 2 of Zll, we 
computed the average radius of the hottest zone, and its standard deviation, for the last 200 s 
and 100 s preceding ignition. For the 4.34 km simulation, over the last 200 s, the average 
hot spot radius and standard deviation are 54.0 km and 22.1 km (as compared to 52.3 km 
and 25.5 km for the 8.68 km simulation). Over the last 100 s, the average hot spot radius 
and standard deviation are 54.7 km and 22.5 km (as compared to 54.7 km and 27.3 km for 
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the 8.68 km simulation). This tells us two things. First, the hot spot statistics do not seem 
to change much whether we consider the final 200 s or 100 s preceding ignition. Second, the 
results are very similar to the 8.68 km simulation, implying that we have sufficient resolution 
in our determination of likely ignition radii. 

Next, as in Zll, we break the final approach to ignition into small time intervals and 
look at properties of the fiow within each time interval. We consider the last 200 s preceding 
ignition, and use time intervals of At^ist = 1.0 s and 0.5 s. Within each time interval, 
we compute the average radius of the hottest cell, the average temperature of the hottest 
cell, and the average radial velocity of the hottest cell. We sort this data into histograms 
to understand the statistics of the last few hot bubbles preceding ignition. In each of the 
following figures, we show histograms for both At^ist = 1.0 s and 0.5 s. Figure 7 shows 
histograms of the hottest cell, sorted by radius, with the colors representing the average 
temperature of the hottest cell over the averaging interval. Figure 8 shows histograms of the 
hottest cell, sorted by radius, with the colors representing the average radial velocity of the 
hottest cell over the averaging interval. Figure 9 shows histograms of the hottest cell, sorted 
by radius, with the colors representing time to ignition. Overall, the results are consistent 
with our observations in Zll, which we now summarize. Some general observations: 

• From each set of histograms, we see that hot spot is most likely to be found between 
40 km and 75 km off-center. This is consistent with both Figure 6 and the histograms 
from Zll. However, we do not see the slight extended tail observed in Zll, which indi- 
cated a slight preference for the hot spots to lie at larger radii within the distribution. 

• For each set of histograms, we observe that the results are essentially the same regard- 
less of whether At^ist = 1.0 s or 0.5 s is used as the averaging interval. 

Some figure-specific observations: 

• In Figure 7, within each temperature interval, the overall shape of the distribution 
appears roughly the same, with a peak slightly greater than 50 km, indicating that hot 
spots of all temperatures can appear at any radius in the distribution. 

• In Figure 8, nearly all of the hot spots have an outward radial velocity. Also, there is 
a tendency for the hot spots at larger radii to be associated with larger values of Vr as 
expected, since the fiow will carry them away from the center. 

• In Figure 9, we see a reasonably symmetric distribution for all cases, indicating that 
the hot spot radius does not depend strongly on time to ignition. 



- 15 - 



Next, we include a new histogram where we examine whether the hottest ceU is increas- 
ing or decreasing in temperature. Figure 10 shows histograms of the hottest ceU, sorted by 
radius, with the colors indicating whether the temperature of the hottest cell is increasing or 
decreasing with time. We observe that when the hottest zone is within 40 km of the center 
it is almost always heating up, and when the hottest zone is outside of 75 km it is almost 
always cooling down. This 40 to 75 km range is consistent with the previous histograms. 
Another conclusion is that it seems highly unlikely that ignition will occur outside of 75 km 
since hot cells beyond that radius are most likely cooling down. This is in contrast to our 
result from Zll, where we claimed that 100 km was an outer limit for ignition radii. 



3.3. Hot Spot Analysis 

We are interested in the likelihood of multiple ignition points, so we now take a closer 
look at the dynamics of the last few hot spots preceding ignition. In the diagnostics we have 
presented, we have only been able to track the hottest cell in the simulation. We do not 
know, for example, if there are other hot spots elsewhere in the star that are almost as hot 
as the hottest zone. It is possible that at the time of ignition, there are one or more cells not 
directly connected to the ignition cell that have almost reached the ignition threshold. In a 
multiple ignition scenario, such cells could also ignite very shortly after the initial ignition. 
Since the white dwarf explodes within a few seconds of ignition, a multiple ignition scenario 
would require another ignition point to develop within the early phases of the explosion for 
it to have any meaningful impact. We wish to examine the temperature field very close to 
the ignition time to get a feel for how likely the multiple ignition scenario is. 

We have previously defined a failed ignition point as a spike in the plot of the peak 
temperature vs. time that does not run away. We begin by examining the temperature 
distribution in the star during three failed ignition points preceding ignition. Figure 11 is 
a zoom-in of Figure 3 for the final minutes preceding ignition for the 4.34 km simulation. 
Three failed ignition points preceding ignition are encapsulated within the green, blue, and 
black dotted lines. We will examine the temperature distribution in each of these time ranges 
to see if there are hot spots elsewhere in the star. 

Figure 12 shows contours of temperature within the green dotted time region from Figure 
11. We note that for all subsequent temperature visualizations, the blue dot represents the 
center of the star, and has a diameter of 4.34 km, corresponding to the resolution of the 
simulation. Also, all visualization frames are spaced at 0.2 s time intervals. The main 
observation is that in the frames where an orange contour exists, indicative of a hot spot, 
that there are no other regions in the star with comparable temperature. This implies that 
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if this hot spot were to run away, there would be only one ignition point. Figure 13 shows 
contours of temperature within the blue dotted time region. We do see that in the frames 
where an orange contour exists, there are other hot spots in different regions of the star. 
But as the hottest spot floats away and cools ofl[, the temperature of the other hot spot 
does not increase. Again, this implies only a single ignition point. Figure 14 shows contours 
of temperature within the black dotted time region. This visualization is more like the 
green dotted time region, in that there are no regions of the star with a peak temperature 
comparable to the main hot spot, implying there would be only one ignition point if this hot 
spot were to run away. 

Next, we perform another simulation, beginning at the point of ignition, in which we 
have disabled burning for all cells where T > 8 x 10^ K. This will give us a picture of the 
dynamics of nearby hot bubbles that did not initially ignite. The idea here is to let the initial 
ignition point float away and see if any other hot spots reach the ignition temperature soon 
afterwards. The peak temperature in this new simulation is given by the black solid line 
in Figure 11. The visualization of the temperature fleld within the pink dotted time region 
from Figure 11 is shown in Figure 15. We see that the hot bubble containing the ignition 
cell floats away from the center of the star and cools off (because the burning is disabled) as 
it breaks up. More importantly, none of the other hot bubbles not connected to the ignition 
cell increase in temperature to the point of ignition. Altogether, our analysis of the last few 
hot spots does not seem to support multiple ignition points, implying that this scenario is 
much less likely than a single ignition point. 

A caveat to this analysis is that our resolution is still several kilometers. It is possible 
that if one could increase the resolution far beyond what is possible today, even with AMR, 
that many smaller hot spots could exist and the dynamics would be different. 

3.4. Convective Flow Pattern 

In Z09, Zll, we provided visualizations of the convective flow fleld, noting the dipole 
feature observed in non-rotating white dwarfs. We recall that the convectively unstable 
region encompasses only an inner fraction (r ^1030 km) of the star. Outside of this, the flow 
is stable against convection and dominated by large-scale structures with high circumferential 
velocities and a smaller radial component. Figure 16 shows visualizations of the 8.68 km, 
4.34 km, and 2.17 km flow flelds in the convective region at t = 10380 s. As before, we show 
contours of outward and inward radial velocity, as well as contours of increasing burning rates. 
As expected, the burning is strongest near the center of the star. Now, we see the effect 
that resolution has on visualization of the velocity contours. Both the large-scale nature of 
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the flow as well as the smaller-scale eddies are much clearer with increasing resolution. This 
allows us to characterize the flow fleld as a plume, with a small solid angle region and strong 
outward velocity, with a lower speed recirculation. This is in contrast to a dipole, where we 
would expect the magnitude of the outflow and inflow to be more similar. In Figure 17, we 
highlight the plume structure by showing the same 2.17 km flow fleld in more detail, where 
each frame represents a 40 degree rotation from the previous. 

In Figure 18 we observe that, for the 4.34 km simulation at the time of ignition, the 
ignition point lies in a region with positive Vr (consistent with our earlier report that Vr = 
9.5 km s~^), and is almost aligned with the strongest outward plume. We expect this hot 
ignition point to accelerate radially outward to a signiflcant fraction of the sound speed 
within a small fraction of a second; this does not give the parcel of fluid at the ignition point 
enough time to change direction and align exactly with the strong outward plume. 

To get an idea of the structure of the flow outside of the convective region, we visualize 
the flow fleld in the x-y plane. In Figure 19, we plot the radial velocity (U • e^), as well 
as one component of the circumferential velocity, U • 651, where Cq is the unit vector in the 
azimuthal direction in the x-y plane. Both plots use the same scale for positive and negative 
velocities, so we see that the circumferential velocities outside the convective region are 
generally larger than the radial velocities within the core. These circumferential velocities 
in the stably stratifled region may become important in explosion simulations in that they 
may deform hot bubbles or flames passing through that region. 

3.5. Turbulence Structure 

Predictive models for SNe la, in particular turbulent flame models, depend critically 
on the structure of the turbulence in the star. In this section, we use the simulations to 
examine this structure and extract estimates for the turbulent intensity and the integral 
length scale. This will help us understand the state of the turbulence that exists at the start 
of the explosion phase. 

The vast majority of literature on turbulence theory deals with flows that are assumed 
to have constant density. In the present context, the signiflcant variation in density due to 
stratiflcation cannot be ignored, and has to be dealt with carefully. Following von Weizsacker 
(1951), Fleck (1983, 1996) advocated casting the energy balance equation in terms of energy 
density (energy per unit volume) as opposed to speciflc energy (energy per unit mass), 
and we note that the difference is inconsequential for constant density turbulence. Thus, 
the fundamental quantity relevant to the inertial range of a turbulent energy cascade is 
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the energy dissipation rate per unit volume Sy^ which should be expected to scale with 
Sy ^ pu^/l^ where u is the turbulent intensity (rms velocity fluctuation), and / is the integral 
length scale. Subsequently, Kritsuk et al. (2007) used numerical simulation of compressible 
turbulence to demonstrate that an appropriately density-weighted velocity spectrum obeys 
a Kolmogorov-type flve-thirds decay law. Consequently, we consider a density-weighted 
velocity fleld, 

V„(a;) = p"U [g" s-i], (10) 

and its Fourier transform, 

= T[\n{x)] [g" cm^-'- s-i]. (11) 
We then deflne a generalized energy density spectrum as 

En{^) = U ^V„(/.) ■ VUk) dS [g2" cm^-e- (12) 

where Q is the volume of the domain in physical space, the domain of the integral, S(/^), is the 
spherical surface deflned by = and * denotes the complex conjugate. This generalized 
energy density spectrum can only be made dimensionless using Sy and ^ for n = 1/3, 
resulting in the dimensionless group £y^^^K.^^^Ei/s{i<i). Therefore, plotting £"1/3 (/^) should 
present a flve-thirds decay. Henceforth, we only present energy density spectra appropriately 
weighted and omit the 1/3 suflix. We also note that only n = 1/2 corresponds to a real energy 
density. 

In the diagnostics in this section, we consider the local velocity, U, rather than the total 
velocity, U + WoCr. In Zll, we showed that the maximum magnitude of Wq at ignition is 
^0.013 km s~^, so the effect of Wq is insigniflcant on the scales we are interested in. 

Energy density spectra from the 2.17 km simulation at t = 10380 s are shown in Figure 
20(a). The density- weighted velocity fleld has been decomposed into different components, 
speciflcally, the Cartesian components, ^4, Vy^ and 14, and spherical polar components 
(using the convention that 9 is the azimuthal angle with respect to the x-y plane and (f is 
the inclination angle measured from the z-axis), 

, . xV^ + yVy + zV^ -yV^ + xVy xzV^ + yzVy - R^V^ 

r R rR 

where = x^ + y^ + z^ and R^ = x^ + y^. Three energy density spectra are plotted: flrst, the 
mean of the Cartesian components (individual components do not differ signiflcantly from 
that shown); second, the radial component; and third, the circumferential component. This 
decomposition demonstrates that there is signiflcantly less energy in the radial component 
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than in the other components. This is due to the large circumferential velocities in the layers 
outside the convection zone; although the density is lower here, the volume is sufficiently 
large that the resulting energy has a significant contribution to the spectrum. It also appears 
that the radial component decays with an exponent close to five-thirds (if slightly smaller), 
and the other components have a slightly higher exponent. 

To circumvent the issue of large circumferential velocities in the stably stratified region, 
and to remove the signal from the coarse-fine interfaces at wavenumbers around 1152 and 
576, the energy density spectra of a subdomain were constructed. This was achieved by 
applying a smoothing function to the velocity and density fields in such a way that the data 
outside the convection zone were set to zero. Specifically, each field was multiplied by the 
hyperbolic tangent function (1 — tanh[(r — ro)/5])/2, where Tq = 875 km and S = 30 km. 
All of the resulting non-trivial data were at the finest AMR level, and the resulting energy 
density spectra are shown in Figure 20(b). Now, each spectrum collapses to a single curve, 
especially for /^^20, which corresponds to a length of about 250 km. The decay exponent 
of each spectrum is close to five-thirds and presents the characteristic "bump" between the 
inertial and dissipation ranges expected from developed homogeneous isotropic turbulence 
(e.g. Saddoughi & Veeravalli 1994; Porter et al. 1994; Kaneda et al. 2003; Aspden et al. 
2008b). 

To explore the effect of resolution on the turbulence in the convective core. Figures 
20(c) and (d) present the total kinetic energy density spectra for the three resolutions, first 
without scaling (c), and then scaled (d). The spectra are scaled according to computational 
cell width and in keeping with a constant energy dissipation rate. Specifically, the 4.34 km 
simulation spectrum is shifted to higher wavenumbers by a factor of 2, and to lower energy 
density by a factor of 2~^/^, and the 8.68 km spectrum has been shifted by factors of 4 and 
^-5/3^ respectively. The unsealed spectra demonstrate that the large scales are independent 
of resolution (as expected), in the sense that increasing the resolution does not lead to an 
increased level of turbulent intensity. This kind of convective motion is not dominated by 
small-scale processes, and integral quantities are well-captured even at moderate resolutions. 
The 8.68 km simulation has a short inertial range, but this is more extensive at higher 
resolutions. The scaled spectra demonstrate that the dissipation range depends on the 
computational cell width as expected from an ILES-type simulation. In particular, there is 
an effective Kolmogorov length scale that is a function of the cell width; the collapse is not 
exact, but is consistent with previous work, see Aspden et al. (2008b), for example, which 
also contains further discussion of the ILES approach and dependence on resolution. 

The rms velocity in the convective core (r < 875 km), was found by direct mea- 
surement to be approximately 14 km s~^ (the data ranged from 12 km s~^ to 18 km s~^ 
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depending on component and resolution); note that no density weighting was used. Even 
though this estimate is smaUer than previously suggested (^100-500 km s~-^), we argue that 
it is actually an upper bound for the turbulence produced by convection because it includes 
large-scale plume-like flow, which artiflcially inflates the estimate. 

To determine the integral length scale in the convective core, the longitudinal correlation 
functions were evaluated for the Cartesian components of the velocity fleld, along with the 
correlation functions of the radial velocity in each Cartesian direction, where the (second- 
order) velocity correlation function (two-point, one-time) is deflned as 



where r denotes the separation vector. The integral length scale in the x direction, for 
example, is then deflned as the integral of the longitudinal velocity correlation function 



The correlation functions were evaluated both for the density- weighted and non- weighted 
velocities (by replacing Ui by Vi in Equation (14) and the appropriate normalization factor 
in Equation (15), and are compared in Figure 21 by solid and dashed lines, respectively. 
The weighted and non- weighted correlation functions are in close agreement, suggesting that 
measuring the integral length scale is not aflFected by the variations in density. By integrat- 
ing each correlation (ignoring the negative parts) , integral length scales for each component 
were evaluated and are shown by the vertical lines with the corresponding line style. The x 
component appears not be consistent with the other components, probably because there is 
a large plume-like structure roughly aligned with the x-axis. Taking this to be an outlier, the 
mean integral length scale was found to be approximately 169 km (with a standard deviation 
of approximately 8.4 km). 

Averages and standard deviations of integral length scale and rms velocity were eval- 
uated using seven time points over 350 s at the 4.34 km resolution, and were found to be 
approximately 200±50 km and 16±3 km s~^, respectively. 

Taking the integral length scale to be 200 km and the turbulent intensity to be 16 km s~^, 
the speciflc energy dissipation rate e = u^/l is approximately 2 x 10^^ cm^ s~^. The cor- 
responding estimates that were suggested to be necessary for a spontaneous detonation by 
Woosley et al. (2011) were 10 km and 500 km s~^, respectively (see also Lisewski et al. 2000; 
Ropke et al. 2007a; Timmes & Woosley 1992). This gives £ ^ 10^^ cm^ s~^, six orders of 
magnitude larger. The present simulations suggest that the turbulent intensity required for 
a spontaneous detonation can not be produced by convection within the core. 
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4. Conclusions and Discussion 

Overall, our high-resolution simulations agreed with the findings of Zll regarding the 
ignition radius of 50 km with a likely range of 40 km to 75 km. We do note that the outer 
limit of 100 km reported in Zll is probably too large, as we do not see any hot bubbles 
at that radius that are still increasing in temperature. By looking closely at the dynamics 
of the last few hot spots, we conclude that the multiple ignition scenario is unlikely. With 
improved resolution, we now describe the large-scale coherent structure in the convective 
field as a plume, rather than a jet, and have a better understanding of the turbulent nature 
of the fiow. 

These findings, together with those from Zll, indicate that a single-point, oflF-center 
ignition is the most likely scenario for SNe la. At the radii we find ignition to be most likely, 
the initial fiame will fioat away faster than it can burn toward the center (see e.g. Plewa et al. 
2004; Zingale & Dursi 2007), making for an asymmetric explosion. This scenario has been 
explored in explosion calculations, potentially giving rise to the "gravitationally confined 
detonation" (Plewa et al. 2004; Jordan et al. 2008), although other groups suggest that this 
mechanism may not be robust (Ropke et al. 2007b). If a single oflF-center ignition fails to blow 
up the star, then it is possible that we would need to wait for the next ignition point, perhaps 
tens of seconds later, or cycle through many widely spaced ignitions until we ignite closer 
to the center (i.e. many successive false starts). Alternately, some type of pulsational model 
may ensue (Ivanova et al. 1974; Khokhlov 1991; Bravo & Garcia-Senz 2006). With these 
results, the challenge to the explosion modelers is to demonstrate that the single-degenerate 
Chandrasekhar mass white dwarf model can produce robust explosions resulting from single- 
point, off-center ignition. Observations may show support for asymmetric models (Maeda 
et al. 2010), but some radiative transfer calculations seem to preclude extreme amounts of 
asymmetry (Blondin et al. 2011). 

We conclude by summarizing the various components of the convecting white dwarf and 
give characteristic length and velocity scales for each; Figure 22 presents this information 
in schematic form. Buoyancy drives a large-scale fiow in the convective core, which extends 
to a radius on the order of 1000 km. This large-scale fiow is composed of plumes around 
100 km wide and several hundred km long with a bulk velocity around 100 km s~^. These 
plumes drive turbulence in the core with an rms velocity and integral length scale that 
were estimated to be on the order of 16 km s~^ and 200 km, respectively. This level of 
turbulence is far below that required for a spontaneous detonation to occur. The stably 
stratified region outside the convective core, extending from ^1000 km to ^1900 km, is 
made up of circumferential shear layers, with a smaller radial velocity component. These 
shear layers are on the order of 100 km deep, several hundred km long, with typical velocities 
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on the order of 100 km and peak velocities that may be in excess of 250 km s"-*^. The 
burning of a single off-center ignition would be dominated at early times by the laminar 
flame speed (on the order of 50 km s~^), and the level of turbulence in the core is unlikely to 
deform the flame very much at all. Furthermore, Aspden et al. (2011) found that large-scale 
entrainment was the dominant process in the evolution of a burning bubble, and that the 
flame speed (turbulent or laminar) even up to 100 km s~^ did not signiflcantly affect the 
evolution. Therefore, the turbulence produced by convection in the core is unlikely to play a 
signiflcant role in the explosion. As the bubble reaches the edge of the convective core, it will 
be ^500 km across moving with a rise speed on the order of 1000 km s~^. The turbulence 
within the bubble itself is likely to have an rms velocity on the order of 100 km s~^ on an 
integral length scale of a few tens of km. In the past, it has been suggested that the convective 
boundary lies at the density suggested for a deflagration-to-detonation transition (Piro & 
Chang 2008). Although the velocities in the core are unlikely to affect the bubble as it rises, 
the circumferential velocities in the stable region are much greater and may interact strongly 
with the bubble as it passes through this region. We plan to investigate this interaction in 
future work. 
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Fig. 1. — (Left) For data on the Cartesian grid (shown here in two dimensions), we use a 
cell-centered convention to denote the average value over the computational cell. (Center) 
The base state variables live on a one-dimensional radial array, and can live at cell centers 
or edges. (Right) A graphical depiction of how the base state and Cartesian grid are related. 
Note that there is no direct alignment between the radial cell centers and the Cartesian grid 
cell centers. 



Fig. 2. — Grid structure for our three-level simulations. The base grid has 576^ grid cells 
(8.68 km resolution), and the refined grids have effective 1152^ (4.34 km) and 2304^ (2.17 km) 
grid cells. The red, green, and blue outlines indicate boxes which can contain up to 64^ grid 
cells. (Left) The shaded region indicates the edge of the star, defined by the location where 
p = 10^ g cm~^ at r 1890 km. (Right) In this zoom-in, the shaded region indicates 
the edge of the convective region, defined by the location where p ^ 1.26 x 10^ g cm~^ at 
r ^ 1030 km. The finest grids contain the entire convective region. 
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Fig. 3. — Peak temperature leading up to ignition for the 8.68 km, 4.34 km, and 2.17 km 
simulations. The inset plot shows the long-time behavior of the 8.68 km simulation originally 
presented in Zll. 
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Fig. 4. — Peak Mach number leading up to ignition for the 8.68 km, 4.34 km, and 2.17 km 
simulations. The inset plot shows the long-time behavior of the 8.68 km simulation originally 
presented in Zll. 
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Fig. 5. — Peak radial velocity leading up to ignition for the 8.68 km, 4.34 km, and 2.17 km 
simulations. The inset plot shows the long-time behavior of the 8.68 km simulation originally 
presented in Zll. 
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Fig. 6. — Radial location of the hottest cell as a function of time for the 4.34 km simulation. 
Only the last 200 s before ignition are shown. Here we see that right up to the end of 
the calculation the hot spot location changes rapidly. The horizontal horizontal dashed line 
indicates the average radial position of the hot spot from 200 s to 1 s before ignition. 
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Fig. 7. — Histograms of the hottest ceU, sorted by radius, with the colors representing the 
average temperature of the hottest ceU over the averaging interval for the 4.34 km simulation 
with (Left) Athist = 1.0 5 and (Right) At^ist = 0.5 5. 
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Fig. 8. — Histograms of the hottest ceU, sorted by radius, with the colors representing 
the average radial velocity of the hottest cell over the averaging interval for the 4.34 km 
simulation with (Left) At^ist = 1.0 5 and (Right) At^ist = 0.5 s. 




Fig. 9. — Histograms of the hottest ceU, sorted by radius, with the colors representing time 
to ignition for the 4.34 km simulation with (Left) At^ist = 1.0 5 and (Right) At^ist = 0.5 s. 
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Fig. 10. — Histograms of the hottest ceU, sorted by radius, with the colors indicating whether 
the temperature of the hottest ceU is increasing or decreasing with time for the 4.34 km 
simulation with (Left) At^ist = 1.0 5 and (Right) At^ist = 0.5 s. 
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Fig. 11. — Peak temperature during the ^200 s preceding ignition for the 4.34 km simulation. 
The dashed vertical lines indicate time ranges where we will examine whether there are 
multiple hot spots. The inset plot is a zoom-in of the final ^5 s preceding ignition. The 
black curve follows the maximum temperature for a simulation where we disable burning in 
aU cells with T > 8 x 10^ K. 
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Fig. 12. — Temperature contours from t = 10404.0 s to t = 10407.8 s (corresponding to 
the green dotted time range in Figure 11) spaced at 0.2 s time intervals. The contours 
are surfaces indicating where T = 7.15 x 10^ K (green), T = 7.2 x 10^ K (yeUow), and 
T = 7.25 X 10^ K (orange). The blue dot is at the center of the star, and has a diameter of 
4.34 km, which corresponds to the grid cell width for this simulation. 
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Fig. 13. — Temperature contours from t = 10454.6 s to t = 10458.4 s (corresponding to the 
blue dotted time range in Figure 11) spaced at 0.2 s time intervals. The contours are surfaces 
indicating where T = 7.24 x 10^ K (green), T = 7.31 x 10^ K (yellow), and T = 7.38 x 10^ K 
(orange). The blue dot is at the center of the star, and has a diameter of 4.34 km, which 
corresponds to the grid cell width for this simulation. 
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Fig. 14. — Temperature contours from t = 10559.2 s to t = 10560.6 s (corresponding to 
the black dotted time range in Figure 11) spaced at 0.2 s time intervals. The contours 
are surfaces indicating where T = 7.48 x 10^ K (green), T = 7.54 x 10^ K (yellow), and 
T = 7.6 X 10^ K (orange). The blue dot is at the center of the star, and has a diameter of 
4.34 km, which corresponds to the grid cell width for this simulation. 
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Fig. 15. — Temperature contours from t = 10562.0 s to t = 10565.0 s (corresponding to the 
pink dotted time range in Figure 11) spaced at 0.2 s time intervals. The contours are surfaces 
indicating where T = 7.5 x 10^ K (green), T = 7.7 x 10^ K (yeUow), and T = 7.9 x 10^ K 
(orange). The blue dot is at the center of the star, and has a diameter of 4.34 km, which 
corresponds to the grid cell width for this simulation. 
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Fig. 16. — Contours of nuclear energy generation rate (yellow to green to purple, correspond- 
ing to 4 X 10^^, 1.27 X 10^^, and 4 x 10^^ erg g"^ s~^) and radial velocity (red is outflow, 
corresponding to 3 x 10^ and 6 x 10^ cm s~^; blue is inflow, corresponding to —3 x 10^ and 
—6 X 10^ cm s~^) for the (clockwise, from top-left) 8.68 km, 4.34 km, and 2.17 km simulations 
at t= 10380 s. Only the inner r = 1000 km are shown. 
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Fig. 17. — Same data as Figure 16, but here we only show the 2.17 km grid ceU simulation, 
and each image represents a view rotation of 40 degrees of the data from t= 10380 s. 
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Fig. 18. — Planar slice from the 4.34 km simulation at the time of ignition, oriented so the 
center of the star (black dot), the ignition location (green dot), and the center of the strongest 
outward plume lie in the plane. The dots each have a radius of 20 km. Red corresponds to 
Vr > 60 km and blue corresponds to Vr < —60 km s~^. Only the inner r = 1000 km is 
shown. 



-41 - 




Fig. 19. — (Top) Plot of radial velocity (U • e^) in the x-y plane from the 2.17 km simulation 
at t = 10380 s. (Bottom) Plot of U • 651 in the x-y plane from the same dataset. In both 
plots, red = +100 km s~^ and blue = -100 km s~^. The outer dark contour indicates the 
edge of the star, where po ^ 1 x 10^ g cm~^ (r 1030 km). The inner dark contour indicates 
the edge of the convective region, where po ^ 1-26 x 10^ g cm~^ (r 1890 km). 
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Fig. 20. — (a) Energy density spectra of the entire domain at 2.17 km resolution for t = 
10380. Note how the radial spectrum is much lower than the other components. This is due 
to the large circumferential velocities outside the convective core, (b) Energy density spectra 
for the convective core. Note how the curves have collapsed to a single profile, especially 
for /^^20, corresponding to about 250 km. (c) Comparison of the energy density spectra 
at the three different resolutions, (d) As (c), but scaled to demonstrate that the effective 
Kolmogorov length scale is proportional to the computational cell width. 
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Fig. 21. — Longitudinal correlation functions for the turbulence in the convective core at 
2.17 km resolution for t = 10380. Density- weighted and non- weighted correlation functions 
are shown by solid and dashed lines, respectively. The x component presents a larger corre- 
lation because there is a plume-like structure roughly aligned with the x-axis. The integral 
length scales are denoted by the vertical lines of the corresponding color. 




Fig. 22. — Cartoon showing the various features with associated velocities and length scales 
in the white dwarf at the end of convection/start of flame propagation. 
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